Workflow: Capital Asset Pricing Model for Ecological Portfolios
Mark D. Scheuerell
Fish Ecology Division, Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Seattle, WA USA, mark.scheuerell@noaa.gov
0.0.1 Version
This is version 0.17.08.01.
1 Initialization
library(readr)
library(reshape2)
library(MARSS)2 Data
For these analyses we will use time series of counts of adult spring/summer Chinook salmon from the Snake River basin in Oregon and Idaho. This group of fish forms an Evolutionarily Significant Unit, which was listed as threatened under the US Endangered Species Act in 1992.
The data were compiled largely through the efforts of Tom Cooney (Northwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, Portland, OR USA, tom.cooney@noaa.gov).
2.1 Assets
The assets data file has the following columns of information:
pop= abbreviated name of the populationcode= 5-letter abbreviation for MPG (2 letters) & population (3 letters)mpg= name of 1 of 5 Major Population Groupsha= estimated area (hectares) of available spawning habitatcategory= hatchery supplementation category (“impact” or “control”)brood_yr= year when spawning occurred
nS= total number of spawning fishphos= proportion of hatchery-origin spawners
rawdata <- read.csv("https://raw.githubusercontent.com/mdscheuerell/CAPM/master/data/srss_chin.csv")
head(rawdata)
## pop code mpg ha category brood_yr nS phos
## 1 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1955 122 0
## 2 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1956 2606 0
## 3 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1957 1754 0
## 4 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1958 486 0
## 5 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1959 712 0
## 6 CathCr GRCAT Grande Ronde / Imnaha 29.1522 impact 1960 3161 0The first thing to do is remove hatchery-origin spawners from the counts, and then convert the counts to log-density based on the estimated spawning area for each population. In a couple of years, incomplete censuses with small sample sizes led to estimated fractions of wild fish equal to zero. I will treat those as NA.
dat <- rawdata[,c("mpg", "code", "pop", "category", "brood_yr")]
## calculate log-density
dat$Sw <- round(log(rawdata$nS*(1-rawdata$phos)/rawdata$ha), 2)
## replace possible -Inf with NA
dat$Sw[is.infinite(dat$Sw)] <- NAFor easier bookkeeping, I’ll also switch the Imnaha’s MPG over to the Grande Ronde.
## change Imnaha name for easier sorting later
dat[,"code"] <- sub("IRMAI", "GR_IR", dat[,"code"], fixed=TRUE)
## MPG abbreviations
mpg_ID <- c("GR", "MF", "SF", "SR")Only a few of the populations have observations going all the way back to the 1950s, so I’ll trim the data accordingly.
## use generation time of 4 years
gt <- 4
## first brood year to consider
yr_first <- 1960
## last brood year to consider
yr_last <- 2011
## years of interest
t_index <- seq(yr_first, yr_last-gt)
## length of time period
TT <- length(t_index)
## subset data
dat <- subset(dat, brood_yr>=yr_first & brood_yr<=yr_last)I will use the MARSS package to fit the CAPM model, which requires the data to be formatted as an [N x T] matrix.
## reshape assets
assets <- dcast(dat, brood_yr ~ code, value.var="Sw")[,-1]
## compute differences
assets <- t(apply(assets, 2, diff, lag=gt))
matplot(t(assets), type="b")
## number of popns
NN <- dim(assets)[1]
## names of popns
names_pop <- rownames(assets)Finally, here are some nicer, longer names for the various populations.
# set up naming & ordering scheme
longnames <- c(
"Wenaha",
"Grande Ronde",
"Catherine",
"Minam",
"Lostine",
"Imnaha",
"SF Salmon",
"Secesh",
"EFSF Salmon",
"Big",
"Sulphur",
"Bear Valley",
"Marsh",
"Loon",
"Camas",
"Yankee Fork",
"Valley",
"UM Salmon",
"LM Salmon",
"EF Salmon",
"Lemhi")
name_ord <- c(6,3,5,4,2,1,12,10,15,14,13,11,9,7,8,20,21,19,18,17,16)
names_mat <- data.frame(ord=name_ord, ID=names_pop)
names_mat <- names_mat[order(names_mat$ord),]
names_mat$long <- longnames2.2 Market index
We used the monthly PDO data provided by the University of Washington’s Joint Institute for the Study of the Atmosphere and Ocean (JISAO), which are available here. We begin by downloading the raw PDO data and viewing the metadata.
## raw PDO data from UW
PDO_raw <- readLines("http://research.jisao.washington.edu/pdo/PDO.latest")
## line with data headers
hdr_pdo <- which(lapply(PDO_raw,grep,pattern="YEAR")==1, arr.ind=TRUE)
## print PDO metadata
print(PDO_raw[seq(hdr_pdo-1)],quote=FALSE)## [1]
## [2] PDO INDEX
## [3]
## [4] If the columns of the table appear without formatting on your browser, use http://research.jisao.washington.edu/pdo/PDO.latest.txt
## [5]
## [6] Updated standardized values for the PDO index, derived as the
## [7] leading PC of monthly SST anomalies in the North Pacific Ocean,
## [8] poleward of 20N. The monthly mean global average SST anomalies
## [9] are removed to separate this pattern of variability from any
## [10] "global warming" signal that may be present in the data.
## [11]
## [12]
## [13] For more details, see:
## [14]
## [15] Zhang, Y., J.M. Wallace, D.S. Battisti, 1997:
## [16] ENSO-like interdecadal variability: 1900-93. J. Climate, 10, 1004-1020.
## [17]
## [18] Mantua, N.J. and S.R. Hare, Y. Zhang, J.M. Wallace, and R.C. Francis,1997:
## [19] A Pacific interdecadal climate oscillation with impacts on salmon
## [20] production. Bulletin of the American Meteorological Society, 78,
## [21] pp. 1069-1079.
## [22] (available via the internet at url:
## [23] http://www.atmos.washington.edu/~mantua/abst.PDO.html)
## [24]
## [25] Data sources for this index are:
## [26] UKMO Historical SST data set for 1900-81;
## [27] Reynold's Optimally Interpolated SST (V1) for January 1982-Dec 2001)
## [28] *** OI SST Version 2 (V2) beginning January 2002 -
## [29]
## [30] missing value flag is -9999
## [31]
## [32]Next, we will extract the actual PDO indices for the years of interest and inspect the file contents.
## PDO data for years of interest
dat_PDO <- read.table("http://research.jisao.washington.edu/pdo/PDO.latest",
header=FALSE, stringsAsFactors=FALSE,
skip=hdr_pdo + yr_first - 1900 + 1 + 2, nrows=TT)
colnames(dat_PDO) <- unlist(strsplit(tolower(PDO_raw[hdr_pdo]), split="\\s+"))
dat_PDO## year jan feb mar apr may jun jul aug sep oct nov dec
## 1 1962 -1.29 -1.15 -1.42 -0.80 -1.22 -1.62 -1.46 -0.48 -1.58 -1.55 -0.37 -0.96
## 2 1963 -0.33 -0.16 -0.54 -0.41 -0.65 -0.88 -1.00 -1.03 0.45 -0.52 -2.08 -1.08
## 3 1964 0.01 -0.21 -0.87 -1.03 -1.91 -0.32 -0.51 -1.03 -0.68 -0.37 -0.80 -1.52
## 4 1965 -1.24 -1.16 0.04 0.62 -0.66 -0.80 -0.47 0.20 0.59 -0.36 -0.59 0.06
## 5 1966 -0.82 -0.03 -1.29 0.06 -0.53 0.16 0.26 -0.35 -0.33 -1.17 -1.15 -0.32
## 6 1967 -0.20 -0.18 -1.20 -0.89 -1.24 -1.16 -0.89 -1.24 -0.72 -0.64 -0.05 -0.40
## 7 1968 -0.95 -0.40 -0.31 -1.03 -0.53 -0.35 0.53 0.19 0.06 -0.34 -0.44 -1.27
## 8 1969 -1.26 -0.95 -0.50 -0.44 -0.20 0.89 0.10 -0.81 -0.66 1.12 0.15 1.38
## 9 1970 0.61 0.43 1.33 0.43 -0.49 0.06 -0.68 -1.63 -1.67 -1.39 -0.80 -0.97
## 10 1971 -1.90 -1.74 -1.68 -1.59 -1.55 -1.55 -2.20 -0.15 0.21 -0.22 -1.25 -1.87
## 11 1972 -1.99 -1.83 -2.09 -1.65 -1.57 -1.87 -0.83 0.25 0.17 0.11 0.57 -0.33
## 12 1973 -0.46 -0.61 -0.50 -0.69 -0.76 -0.97 -0.57 -1.14 -0.51 -0.87 -1.81 -0.76
## 13 1974 -1.22 -1.65 -0.90 -0.52 -0.28 -0.31 -0.08 0.27 0.44 -0.10 0.43 -0.12
## 14 1975 -0.84 -0.71 -0.51 -1.30 -1.02 -1.16 -0.40 -1.07 -1.23 -1.29 -2.08 -1.61
## 15 1976 -1.14 -1.85 -0.96 -0.89 -0.68 -0.67 0.61 1.28 0.82 1.11 1.25 1.22
## 16 1977 1.65 1.11 0.72 0.30 0.31 0.42 0.19 0.64 -0.55 -0.61 -0.72 -0.69
## 17 1978 0.34 1.45 1.34 1.29 0.90 0.15 -1.24 -0.56 -0.44 0.10 -0.07 -0.43
## 18 1979 -0.58 -1.33 0.30 0.89 1.09 0.17 0.84 0.52 1.00 1.06 0.48 -0.42
## 19 1980 -0.11 1.32 1.09 1.49 1.20 -0.22 0.23 0.51 0.10 1.35 0.37 -0.10
## 20 1981 0.59 1.46 0.99 1.45 1.75 1.69 0.84 0.18 0.42 0.18 0.80 0.67
## 21 1982 0.34 0.20 0.19 -0.19 -0.58 -0.78 0.58 0.39 0.84 0.37 -0.25 0.26
## 22 1983 0.56 1.14 2.11 1.87 1.80 2.36 3.51 1.85 0.91 0.96 1.02 1.69
## 23 1984 1.50 1.21 1.77 1.52 1.30 0.18 -0.18 -0.03 0.67 0.58 0.71 0.82
## 24 1985 1.27 0.94 0.57 0.19 0.00 0.18 1.07 0.81 0.44 0.29 -0.75 0.38
## 25 1986 1.12 1.61 2.18 1.55 1.16 0.89 1.38 0.22 0.22 1.00 1.77 1.77
## 26 1987 1.88 1.75 2.10 2.16 1.85 0.73 2.01 2.83 2.44 1.36 1.47 1.27
## 27 1988 0.93 1.24 1.42 0.94 1.20 0.74 0.64 0.19 -0.37 -0.10 -0.02 -0.43
## 28 1989 -0.95 -1.02 -0.83 -0.32 0.47 0.36 0.83 0.09 0.05 -0.12 -0.50 -0.21
## 29 1990 -0.30 -0.65 -0.62 0.27 0.44 0.44 0.27 0.11 0.38 -0.69 -1.69 -2.23
## 30 1991 -2.02 -1.19 -0.74 -1.01 -0.51 -1.47 -0.10 0.36 0.65 0.49 0.42 0.09
## 31 1992 0.05 0.31 0.67 0.75 1.54 1.26 1.90 1.44 0.83 0.93 0.93 0.53
## 32 1993 0.05 0.19 0.76 1.21 2.13 2.34 2.35 2.69 1.56 1.41 1.24 1.07
## 33 1994 1.21 0.59 0.80 1.05 1.23 0.46 0.06 -0.79 -1.36 -1.32 -1.96 -1.79
## 34 1995 -0.49 0.46 0.75 0.83 1.46 1.27 1.71 0.21 1.16 0.47 -0.28 0.16
## 35 1996 0.59 0.75 1.01 1.46 2.18 1.10 0.77 -0.14 0.24 -0.33 0.09 -0.03
## 36 1997 0.23 0.28 0.65 1.05 1.83 2.76 2.35 2.79 2.19 1.61 1.12 0.67
## 37 1998 0.83 1.56 2.01 1.27 0.70 0.40 -0.04 -0.22 -1.21 -1.39 -0.52 -0.44
## 38 1999 -0.32 -0.66 -0.33 -0.41 -0.68 -1.30 -0.66 -0.96 -1.53 -2.23 -2.05 -1.63
## 39 2000 -2.00 -0.83 0.29 0.35 -0.05 -0.44 -0.66 -1.19 -1.24 -1.30 -0.53 0.52
## 40 2001 0.60 0.29 0.45 -0.31 -0.30 -0.47 -1.31 -0.77 -1.37 -1.37 -1.26 -0.93
## 41 2002** 0.27 -0.64 -0.43 -0.32 -0.63 -0.35 -0.31 0.60 0.43 0.42 1.51 2.10
## 42 2003** 2.09 1.75 1.51 1.18 0.89 0.68 0.96 0.88 0.01 0.83 0.52 0.33
## 43 2004** 0.43 0.48 0.61 0.57 0.88 0.04 0.44 0.85 0.75 -0.11 -0.63 -0.17
## 44 2005** 0.44 0.81 1.36 1.03 1.86 1.17 0.66 0.25 -0.46 -1.32 -1.50 0.20
## 45 2006** 1.03 0.66 0.05 0.40 0.48 1.04 0.35 -0.65 -0.94 -0.05 -0.22 0.14
## 46 2007** 0.01 0.04 -0.36 0.16 -0.10 0.09 0.78 0.50 -0.36 -1.45 -1.08 -0.58
## 47 2008** -1.00 -0.77 -0.71 -1.52 -1.37 -1.34 -1.67 -1.70 -1.55 -1.76 -1.25 -0.87
## 48 2009** -1.40 -1.55 -1.59 -1.65 -0.88 -0.31 -0.53 0.09 0.52 0.27 -0.40 0.08Notice that values in the year column in the PDO data file contains an ** superscript from 2002 onward to indicate the switch to version 2 of the sea surface temperature (SST) data. We will fix those as well.
## fix 'year' col
dat_PDO$year <- seq(1:TT) + yr_first - 1
## mean PDO from Apr-Oct indexed by brood year
PDO_ts <- apply(dat_PDO[,-1],1,mean)
dat_PDO <- data.frame(year=dat_PDO$year, PDO=PDO_ts)
## market index
market <- round(dat_PDO$PDO, 3)
plot.ts(market)
2.3 Risk-free index
For these purposes, the risk free index is zero for all years. The idea is that any population with a standardized return of zero is replacing itself and therefore not declining.
free <- rep(0, TT)3 CAPM in MARSS
3.1 v1: a static; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## equal var-cov
# QQ[bb,bb] <- "beta.cov"
# diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1301 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1301 iterations.
## Log-likelihood: -372.3136
## AIC: 812.6272 AICc: 822.3414
##
## Estimate
## R.diag 0.5000
## Q.beta1 1.9085
## Q.beta12 1.7668
## Q.beta13 2.5422
## Q.beta14 2.1624
## Q.beta15 2.2498
## Q.beta16 2.3169
## Q.beta2 1.8029
## Q.beta23 2.2418
## Q.beta24 1.8395
## Q.beta25 1.9753
## Q.beta26 2.0533
## Q.beta3 3.4610
## Q.beta34 2.9889
## Q.beta35 3.0685
## Q.beta36 3.1475
## Q.beta4 2.6078
## Q.beta45 2.6534
## Q.beta46 2.7143
## Q.beta5 2.7211
## Q.beta56 2.7902
## Q.beta6 2.8631
## x0.X1 0.0506
## x0.X2 0.1166
## x0.X3 0.2123
## x0.X4 0.1170
## x0.X5 0.0219
## x0.X6 0.0690
## x0.X7 0.2993
## x0.X8 1.5459
## x0.X9 -0.0543
## x0.X10 -0.2095
## x0.X11 -0.3033
## x0.X12 0.1886
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2224 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2224 iterations.
## Log-likelihood: -397.1267
## AIC: 862.2534 AICc: 872.4243
##
## Estimate
## R.diag 0.4909
## Q.beta1 12.9265
## Q.beta12 14.1924
## Q.beta13 11.1113
## Q.beta14 11.7245
## Q.beta15 12.3104
## Q.beta16 11.0612
## Q.beta2 15.6232
## Q.beta23 12.3696
## Q.beta24 13.0052
## Q.beta25 13.5118
## Q.beta26 12.0437
## Q.beta3 10.3303
## Q.beta34 10.6704
## Q.beta35 10.5630
## Q.beta36 9.0572
## Q.beta4 11.0871
## Q.beta45 11.1514
## Q.beta46 9.6881
## Q.beta5 11.7242
## Q.beta56 10.5450
## Q.beta6 9.7275
## x0.X1 0.0847
## x0.X2 0.0837
## x0.X3 -0.1976
## x0.X4 -0.0232
## x0.X5 0.0689
## x0.X6 0.0163
## x0.X7 -0.0795
## x0.X8 0.2604
## x0.X9 -0.3900
## x0.X10 0.0398
## x0.X11 -0.1894
## x0.X12 -0.1460
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 655 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 655 iterations.
## Log-likelihood: -145.5705
## AIC: 317.1409 AICc: 320.1245
##
## Estimate
## R.diag 0.26835
## Q.beta1 1.24348
## Q.beta12 0.98983
## Q.beta13 0.91831
## Q.beta2 0.78821
## Q.beta23 0.72955
## Q.beta3 0.68857
## x0.X1 -0.00702
## x0.X2 -0.07964
## x0.X3 0.13335
## x0.X4 0.67648
## x0.X5 0.75019
## x0.X6 0.93047
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1560 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1560 iterations.
## Log-likelihood: -372.5814
## AIC: 813.1629 AICc: 822.9571
##
## Estimate
## R.diag 0.4367
## Q.beta1 8.7974
## Q.beta12 7.5370
## Q.beta13 6.6868
## Q.beta14 8.8338
## Q.beta15 9.2384
## Q.beta16 11.6973
## Q.beta2 6.4641
## Q.beta23 5.7336
## Q.beta24 7.5611
## Q.beta25 7.8883
## Q.beta26 10.0038
## Q.beta3 5.0860
## Q.beta34 6.7095
## Q.beta35 7.0034
## Q.beta36 8.8786
## Q.beta4 8.8781
## Q.beta45 9.3050
## Q.beta46 11.7648
## Q.beta5 9.8061
## Q.beta56 12.3538
## Q.beta6 15.6003
## x0.X1 0.2527
## x0.X2 0.1038
## x0.X3 0.0971
## x0.X4 0.2392
## x0.X5 0.1699
## x0.X6 0.2439
## x0.X7 0.2062
## x0.X8 0.2658
## x0.X9 0.4366
## x0.X10 -0.1898
## x0.X11 -0.3480
## x0.X12 -0.4036
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.1.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.2 v1: a static; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1380 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1380 iterations.
## Log-likelihood: -378.8426
## AIC: 787.6852 AICc: 789.5034
##
## Estimate
## R.diag 0.6063
## Q.beta.var 2.0867
## Q.beta.cov 2.0866
## x0.X1 0.0302
## x0.X2 0.0503
## x0.X3 0.1612
## x0.X4 0.1015
## x0.X5 -0.0105
## x0.X6 0.0184
## x0.X7 0.1434
## x0.X8 0.3911
## x0.X9 0.2313
## x0.X10 0.5211
## x0.X11 0.0687
## x0.X12 0.4557
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1390 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1390 iterations.
## Log-likelihood: -404.0779
## AIC: 838.1558 AICc: 840.053
##
## Estimate
## R.diag 0.6716
## Q.beta.var 9.1373
## Q.beta.cov 9.1372
## x0.X1 0.0181
## x0.X2 0.0015
## x0.X3 -0.1194
## x0.X4 -0.0641
## x0.X5 0.0122
## x0.X6 -0.0411
## x0.X7 -0.0973
## x0.X8 0.0320
## x0.X9 -0.2713
## x0.X10 -0.1796
## x0.X11 -0.2088
## x0.X12 -0.0799
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1029 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1029 iterations.
## Log-likelihood: -147.174
## AIC: 312.348 AICc: 313.7765
##
## Estimate
## R.diag 0.303214
## Q.beta.var 0.753774
## Q.beta.cov 0.753725
## x0.X1 0.000199
## x0.X2 -0.035161
## x0.X3 0.082217
## x0.X4 0.634399
## x0.X5 0.898669
## x0.X6 0.885139
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1394 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1394 iterations.
## Log-likelihood: -387.9918
## AIC: 805.9837 AICc: 807.8157
##
## Estimate
## R.diag 0.5493
## Q.beta.var 7.4219
## Q.beta.cov 7.4218
## x0.X1 0.2113
## x0.X2 0.0932
## x0.X3 0.1140
## x0.X4 0.2188
## x0.X5 0.1759
## x0.X6 0.0990
## x0.X7 0.0423
## x0.X8 -0.0261
## x0.X9 0.2312
## x0.X10 -0.1561
## x0.X11 0.1660
## x0.X12 -0.2863
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.2.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.3 v1: a static; b diag & unequal
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## lower R block for beta
bb <- aa + nn
## diag & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1276 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1276 iterations.
## Log-likelihood: -420.0602
## AIC: 878.1204 AICc: 881.0434
##
## Estimate
## R.diag 1.014417
## Q.beta1 0.016075
## Q.beta2 0.300899
## Q.beta3 0.096780
## Q.beta4 0.000323
## Q.beta5 0.028867
## Q.beta6 0.012729
## x0.X1 -0.058006
## x0.X2 0.038213
## x0.X3 -0.033504
## x0.X4 -0.018423
## x0.X5 -0.070038
## x0.X6 -0.054581
## x0.X7 -0.080212
## x0.X8 1.108690
## x0.X9 -0.737083
## x0.X10 -0.014390
## x0.X11 -0.107951
## x0.X12 0.222093
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2156 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2156 iterations.
## Log-likelihood: -459.2719
## AIC: 956.5437 AICc: 959.5959
##
## Estimate
## R.diag 1.676423
## Q.beta1 0.000370
## Q.beta2 0.054952
## Q.beta3 0.000115
## Q.beta4 0.012216
## Q.beta5 0.017199
## Q.beta6 0.050396
## x0.X1 -0.086090
## x0.X2 -0.039976
## x0.X3 -0.233259
## x0.X4 -0.059126
## x0.X5 -0.150412
## x0.X6 -0.322697
## x0.X7 -0.192926
## x0.X8 0.684116
## x0.X9 -0.455433
## x0.X10 0.079412
## x0.X11 -0.152422
## x0.X12 -0.205091
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 108 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 108 iterations.
## Log-likelihood: -172.9275
## AIC: 365.8549 AICc: 367.6149
##
## Estimate
## R.diag 0.4606
## Q.beta1 0.4191
## Q.beta2 0.2020
## Q.beta3 0.1118
## x0.X1 -0.0528
## x0.X2 -0.0472
## x0.X3 0.1730
## x0.X4 0.3972
## x0.X5 1.0640
## x0.X6 0.7844
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2756 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2756 iterations.
## Log-likelihood: -479.2527
## AIC: 996.5053 AICc: 999.4511
##
## Estimate
## R.diag 1.64e+00
## Q.beta1 2.11e-04
## Q.beta2 9.32e-05
## Q.beta3 1.01e-04
## Q.beta4 5.99e-05
## Q.beta5 1.37e-01
## Q.beta6 9.92e-01
## x0.X1 -1.09e-01
## x0.X2 -2.27e-01
## x0.X3 -1.74e-01
## x0.X4 -1.01e-01
## x0.X5 -2.90e-01
## x0.X6 -3.54e-01
## x0.X7 -2.20e-01
## x0.X8 -2.92e-01
## x0.X9 -1.36e-02
## x0.X10 -4.21e-01
## x0.X11 -3.34e-01
## x0.X12 -1.21e+00
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.3.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.4 v2: a unconstrained; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j,k] <- paste0("alpha",j,k)
QQ[k,j] <- paste0("alpha",j,k)
}
}
## lower R block for beta
bb <- aa + nn
## equal var-cov
# QQ[bb,bb] <- "beta.cov"
# diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1886 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1886 iterations.
## Log-likelihood: -350.0825
## AIC: 810.165 AICc: 837.665
##
## Estimate
## R.diag 0.21763
## Q.alpha1 0.35201
## Q.alpha12 0.28274
## Q.alpha13 0.46744
## Q.alpha14 0.46376
## Q.alpha15 0.33098
## Q.alpha16 0.34731
## Q.alpha2 0.50670
## Q.alpha23 0.40905
## Q.alpha24 0.25186
## Q.alpha25 0.52711
## Q.alpha26 0.32887
## Q.alpha3 0.64873
## Q.alpha34 0.60303
## Q.alpha35 0.49702
## Q.alpha36 0.44837
## Q.alpha4 0.66726
## Q.alpha45 0.31397
## Q.alpha46 0.42392
## Q.alpha5 0.61457
## Q.alpha56 0.38241
## Q.alpha6 0.39532
## Q.beta1 0.08182
## Q.beta12 0.13085
## Q.beta13 0.06558
## Q.beta14 0.03405
## Q.beta15 -0.00181
## Q.beta16 0.04432
## Q.beta2 0.28132
## Q.beta23 0.05297
## Q.beta24 0.03642
## Q.beta25 0.13262
## Q.beta26 0.07967
## Q.beta3 0.09699
## Q.beta34 0.07430
## Q.beta35 0.01698
## Q.beta36 0.07638
## Q.beta4 0.18435
## Q.beta45 0.53115
## Q.beta46 0.24618
## Q.beta5 2.18821
## Q.beta56 0.80106
## Q.beta6 0.34428
## x0.X1 0.56813
## x0.X2 1.97133
## x0.X3 0.82577
## x0.X4 0.25570
## x0.X5 1.44091
## x0.X6 0.64654
## x0.X7 0.91927
## x0.X8 3.62940
## x0.X9 -0.01827
## x0.X10 -0.08166
## x0.X11 1.00034
## x0.X12 0.38951
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1905 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1905 iterations.
## Log-likelihood: -344.0416
## AIC: 798.0832 AICc: 827.0034
##
## Estimate
## R.diag 0.242256
## Q.alpha1 1.408860
## Q.alpha12 1.613053
## Q.alpha13 0.734787
## Q.alpha14 0.674628
## Q.alpha15 1.414349
## Q.alpha16 1.470110
## Q.alpha2 1.874476
## Q.alpha23 0.882608
## Q.alpha24 0.801018
## Q.alpha25 1.615979
## Q.alpha26 1.648587
## Q.alpha3 0.936178
## Q.alpha34 1.026439
## Q.alpha35 0.719197
## Q.alpha36 0.226078
## Q.alpha4 1.165429
## Q.alpha45 0.656512
## Q.alpha46 0.039200
## Q.alpha5 1.420653
## Q.alpha56 1.493419
## Q.alpha6 2.064064
## Q.beta1 0.001881
## Q.beta12 -0.000151
## Q.beta13 -0.001419
## Q.beta14 -0.013128
## Q.beta15 -0.001147
## Q.beta16 0.022776
## Q.beta2 0.002387
## Q.beta23 -0.003008
## Q.beta24 0.005799
## Q.beta25 0.000618
## Q.beta26 -0.006824
## Q.beta3 0.005940
## Q.beta34 0.004748
## Q.beta35 0.000453
## Q.beta36 -0.012415
## Q.beta4 0.103576
## Q.beta45 0.009616
## Q.beta46 -0.173020
## Q.beta5 0.000974
## Q.beta56 -0.015933
## Q.beta6 0.293174
## x0.X1 -0.391740
## x0.X2 -0.870161
## x0.X3 -0.664801
## x0.X4 -0.190172
## x0.X5 -0.380069
## x0.X6 -0.301613
## x0.X7 -0.530841
## x0.X8 -0.243570
## x0.X9 -0.757735
## x0.X10 -0.268795
## x0.X11 -0.692295
## x0.X12 -0.375913
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1303 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1303 iterations.
## Log-likelihood: -131.9363
## AIC: 301.8725 AICc: 308.4242
##
## Estimate
## R.diag 0.121
## Q.alpha1 0.443
## Q.alpha12 0.357
## Q.alpha13 0.386
## Q.alpha2 0.293
## Q.alpha23 0.289
## Q.alpha3 0.428
## Q.beta1 0.192
## Q.beta12 0.153
## Q.beta13 0.142
## Q.beta2 0.122
## Q.beta23 0.113
## Q.beta3 0.105
## x0.X1 -0.582
## x0.X2 -0.581
## x0.X3 -0.246
## x0.X4 0.243
## x0.X5 0.431
## x0.X6 0.357
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1891 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1891 iterations.
## Log-likelihood: -332.1111
## AIC: 774.2221 AICc: 801.9699
##
## Estimate
## R.diag 0.24632
## Q.alpha1 1.41757
## Q.alpha12 0.85683
## Q.alpha13 1.02300
## Q.alpha14 1.13028
## Q.alpha15 1.19826
## Q.alpha16 1.75236
## Q.alpha2 0.73748
## Q.alpha23 0.66275
## Q.alpha24 0.71257
## Q.alpha25 0.60989
## Q.alpha26 0.89624
## Q.alpha3 0.74942
## Q.alpha34 0.82637
## Q.alpha35 0.84293
## Q.alpha36 1.25107
## Q.alpha4 0.91557
## Q.alpha45 0.94303
## Q.alpha46 1.41794
## Q.alpha5 1.07329
## Q.alpha56 1.57804
## Q.alpha6 2.46076
## Q.beta1 0.02213
## Q.beta12 -0.01139
## Q.beta13 -0.00460
## Q.beta14 0.01833
## Q.beta15 0.08796
## Q.beta16 0.04825
## Q.beta2 0.00599
## Q.beta23 0.00245
## Q.beta24 -0.00945
## Q.beta25 -0.04558
## Q.beta26 -0.02495
## Q.beta3 0.00104
## Q.beta34 -0.00382
## Q.beta35 -0.01853
## Q.beta36 -0.01011
## Q.beta4 0.01521
## Q.beta45 0.07293
## Q.beta46 0.04002
## Q.beta5 0.35056
## Q.beta56 0.19215
## Q.beta6 0.10543
## x0.X1 -0.67205
## x0.X2 -1.00755
## x0.X3 -0.63267
## x0.X4 -0.55048
## x0.X5 -0.43898
## x0.X6 -0.33243
## x0.X7 -0.57318
## x0.X8 -0.52855
## x0.X9 -0.24796
## x0.X10 -0.79321
## x0.X11 -0.57759
## x0.X12 -1.23466
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.4.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.5 v2: a unconstrained; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j,k] <- paste0("alpha",j,k)
QQ[k,j] <- paste0("alpha",j,k)
}
}
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## diagonal & unequal
# diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
# for(j in 1:(nn-1)) {
# for(k in (j+1):nn) {
# QQ[j+nn,k+nn] <- paste0("beta",j,k)
# QQ[k+nn,j+nn] <- paste0("beta",j,k)
# }
# }
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 2003 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2003 iterations.
## Log-likelihood: -358.099
## AIC: 788.198 AICc: 799.1609
##
## Estimate
## R.diag 0.2833
## Q.alpha1 0.3294
## Q.alpha12 0.2799
## Q.alpha13 0.4307
## Q.alpha14 0.4216
## Q.alpha15 0.2948
## Q.alpha16 0.3199
## Q.alpha2 0.6465
## Q.alpha23 0.3446
## Q.alpha24 0.1407
## Q.alpha25 0.4108
## Q.alpha26 0.3044
## Q.alpha3 0.6389
## Q.alpha34 0.6086
## Q.alpha35 0.6259
## Q.alpha36 0.5225
## Q.alpha4 0.6851
## Q.alpha45 0.4276
## Q.alpha46 0.4442
## Q.alpha5 1.3891
## Q.alpha56 0.8241
## Q.alpha6 0.5904
## Q.beta.var 0.0675
## Q.beta.cov 0.0675
## x0.X1 0.8072
## x0.X2 0.0792
## x0.X3 1.5557
## x0.X4 1.6881
## x0.X5 1.0960
## x0.X6 0.8313
## x0.X7 0.9286
## x0.X8 1.6950
## x0.X9 1.0510
## x0.X10 1.1087
## x0.X11 0.6499
## x0.X12 1.0357
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 2053 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2053 iterations.
## Log-likelihood: -348.8623
## AIC: 769.7246 AICc: 781.2073
##
## Estimate
## R.diag 2.90e-01
## Q.alpha1 1.41e+00
## Q.alpha12 1.56e+00
## Q.alpha13 6.73e-01
## Q.alpha14 6.51e-01
## Q.alpha15 1.43e+00
## Q.alpha16 1.54e+00
## Q.alpha2 1.77e+00
## Q.alpha23 8.30e-01
## Q.alpha24 8.10e-01
## Q.alpha25 1.59e+00
## Q.alpha26 1.63e+00
## Q.alpha3 9.73e-01
## Q.alpha34 1.10e+00
## Q.alpha35 6.56e-01
## Q.alpha36 9.85e-02
## Q.alpha4 1.26e+00
## Q.alpha45 6.27e-01
## Q.alpha46 -5.57e-02
## Q.alpha5 1.46e+00
## Q.alpha56 1.60e+00
## Q.alpha6 2.31e+00
## Q.beta.var 7.07e-05
## Q.beta.cov 3.70e-05
## x0.X1 -5.90e-01
## x0.X2 -9.56e-01
## x0.X3 -2.72e-01
## x0.X4 7.24e-02
## x0.X5 -6.33e-01
## x0.X6 -8.93e-01
## x0.X7 -6.42e-01
## x0.X8 -4.56e-01
## x0.X9 -3.37e-01
## x0.X10 -8.24e-02
## x0.X11 -7.66e-01
## x0.X12 -9.58e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1120 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1120 iterations.
## Log-likelihood: -132.3308
## AIC: 294.6615 AICc: 298.6615
##
## Estimate
## R.diag 0.1204
## Q.alpha1 0.6003
## Q.alpha12 0.4713
## Q.alpha13 0.4842
## Q.alpha2 0.3765
## Q.alpha23 0.3567
## Q.alpha3 0.4764
## Q.beta.var 0.0537
## Q.beta.cov 0.0537
## x0.X1 -0.3677
## x0.X2 -0.4026
## x0.X3 0.0302
## x0.X4 0.3698
## x0.X5 0.6479
## x0.X6 0.6231
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 2106 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 2106 iterations.
## Log-likelihood: -339.5016
## AIC: 751.0032 AICc: 762.0571
##
## Estimate
## R.diag 3.11e-01
## Q.alpha1 1.45e+00
## Q.alpha12 8.24e-01
## Q.alpha13 1.01e+00
## Q.alpha14 1.19e+00
## Q.alpha15 1.35e+00
## Q.alpha16 1.85e+00
## Q.alpha2 6.55e-01
## Q.alpha23 6.21e-01
## Q.alpha24 6.73e-01
## Q.alpha25 6.47e-01
## Q.alpha26 8.95e-01
## Q.alpha3 7.17e-01
## Q.alpha34 8.32e-01
## Q.alpha35 9.11e-01
## Q.alpha36 1.26e+00
## Q.alpha4 9.85e-01
## Q.alpha45 1.11e+00
## Q.alpha46 1.55e+00
## Q.alpha5 1.33e+00
## Q.alpha56 1.84e+00
## Q.alpha6 2.59e+00
## Q.beta.var 5.28e-05
## Q.beta.cov 2.85e-05
## x0.X1 -5.14e-01
## x0.X2 -1.14e+00
## x0.X3 -6.34e-01
## x0.X4 -3.91e-01
## x0.X5 -7.59e-02
## x0.X6 -4.51e-02
## x0.X7 -4.54e-01
## x0.X8 -6.30e-01
## x0.X9 -2.76e-01
## x0.X10 -6.67e-01
## x0.X11 -2.45e-01
## x0.X12 -9.47e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.5.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.6 v3: a equalvarcov; b unconstrained
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
QQ[aa,aa] <- "alpha.cov"
diag(QQ[aa,aa]) <- "alpha.var"
## lower R block for beta
bb <- aa + nn
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
for(j in 1:(nn-1)) {
for(k in (j+1):nn) {
QQ[j+nn,k+nn] <- paste0("beta",j,k)
QQ[k+nn,j+nn] <- paste0("beta",j,k)
}
}
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1950 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1950 iterations.
## Log-likelihood: -360.5338
## AIC: 793.0676 AICc: 804.0305
##
## Estimate
## R.diag 0.4892
## Q.alpha.var 0.3616
## Q.alpha.cov 0.3616
## Q.beta1 0.0382
## Q.beta12 0.0574
## Q.beta13 0.0320
## Q.beta14 0.0248
## Q.beta15 0.0367
## Q.beta16 0.0306
## Q.beta2 0.2277
## Q.beta23 -0.0424
## Q.beta24 -0.1070
## Q.beta25 -0.0641
## Q.beta26 -0.0480
## Q.beta3 0.0845
## Q.beta34 0.1130
## Q.beta35 0.1068
## Q.beta36 0.0857
## Q.beta4 0.1634
## Q.beta45 0.1454
## Q.beta46 0.1158
## Q.beta5 0.1356
## Q.beta56 0.1086
## Q.beta6 0.0871
## x0.X1 0.3421
## x0.X2 0.4242
## x0.X3 0.4651
## x0.X4 0.3915
## x0.X5 0.3153
## x0.X6 0.3377
## x0.X7 0.5544
## x0.X8 1.8386
## x0.X9 0.1294
## x0.X10 0.0441
## x0.X11 -0.0581
## x0.X12 0.2787
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1923 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1923 iterations.
## Log-likelihood: -368.2447
## AIC: 808.4893 AICc: 819.9721
##
## Estimate
## R.diag 0.4578
## Q.alpha.var 1.0297
## Q.alpha.cov 1.0297
## Q.beta1 0.3012
## Q.beta12 0.2121
## Q.beta13 -0.2561
## Q.beta14 -0.1346
## Q.beta15 0.2690
## Q.beta16 0.5721
## Q.beta2 0.1590
## Q.beta23 -0.1796
## Q.beta24 -0.0879
## Q.beta25 0.1911
## Q.beta26 0.3956
## Q.beta3 0.2181
## Q.beta34 0.1151
## Q.beta35 -0.2286
## Q.beta36 -0.4870
## Q.beta4 0.0651
## Q.beta45 -0.1190
## Q.beta46 -0.2607
## Q.beta5 0.2406
## Q.beta56 0.5097
## Q.beta6 1.0917
## x0.X1 -0.1585
## x0.X2 -0.1619
## x0.X3 -0.4460
## x0.X4 -0.2845
## x0.X5 -0.1714
## x0.X6 -0.1760
## x0.X7 -0.2652
## x0.X8 0.2122
## x0.X9 -0.6678
## x0.X10 -0.1822
## x0.X11 -0.3212
## x0.X12 -0.3869
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1219 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1219 iterations.
## Log-likelihood: -134.9975
## AIC: 299.9951 AICc: 303.9951
##
## Estimate
## R.diag 0.1664
## Q.alpha.var 0.2850
## Q.alpha.cov 0.2850
## Q.beta1 0.3401
## Q.beta12 0.2731
## Q.beta13 0.1665
## Q.beta2 0.2209
## Q.beta23 0.1282
## Q.beta3 0.0995
## x0.X1 -0.5229
## x0.X2 -0.5911
## x0.X3 -0.3300
## x0.X4 0.2691
## x0.X5 0.3843
## x0.X6 0.4081
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1934 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1934 iterations.
## Log-likelihood: -347.5903
## AIC: 767.1806 AICc: 778.2345
##
## Estimate
## R.diag 0.3472
## Q.alpha.var 0.6749
## Q.alpha.cov 0.6748
## Q.beta1 0.3710
## Q.beta12 0.1475
## Q.beta13 0.2785
## Q.beta14 0.4376
## Q.beta15 0.7304
## Q.beta16 1.1119
## Q.beta2 0.0598
## Q.beta23 0.1144
## Q.beta24 0.1747
## Q.beta25 0.2846
## Q.beta26 0.4466
## Q.beta3 0.2218
## Q.beta34 0.3309
## Q.beta35 0.5285
## Q.beta36 0.8501
## Q.beta4 0.5167
## Q.beta45 0.8579
## Q.beta46 1.3146
## Q.beta5 1.4689
## Q.beta56 2.1651
## Q.beta6 3.3515
## x0.X1 -1.0269
## x0.X2 -1.1166
## x0.X3 -1.0567
## x0.X4 -1.0128
## x0.X5 -1.1300
## x0.X6 -1.1249
## x0.X7 -0.8751
## x0.X8 -0.6575
## x0.X9 -0.4932
## x0.X10 -1.1376
## x0.X11 -1.2226
## x0.X12 -1.9523
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.6.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.7 v3: a equalvarcov; b equalvarcov
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
QQ[aa,aa] <- "alpha.cov"
diag(QQ[aa,aa]) <- "alpha.var"
## lower R block for beta
bb <- aa + nn
## equal var-cov
QQ[bb,bb] <- "beta.cov"
diag(QQ[bb,bb]) <- "beta.var"
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1916 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1916 iterations.
## Log-likelihood: -366.8833
## AIC: 767.7665 AICc: 770.1024
##
## Estimate
## R.diag 0.5887
## Q.alpha.var 0.3221
## Q.alpha.cov 0.3221
## Q.beta.var 0.0502
## Q.beta.cov 0.0501
## x0.X1 0.4183
## x0.X2 0.4331
## x0.X3 0.5491
## x0.X4 0.4894
## x0.X5 0.3806
## x0.X6 0.4096
## x0.X7 0.4309
## x0.X8 0.6753
## x0.X9 0.5201
## x0.X10 0.8094
## x0.X11 0.3756
## x0.X12 0.7383
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 1937 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1937 iterations.
## Log-likelihood: -376.9972
## AIC: 787.9943 AICc: 790.4326
##
## Estimate
## R.diag 6.26e-01
## Q.alpha.var 1.07e+00
## Q.alpha.cov 1.07e+00
## Q.beta.var 1.10e-04
## Q.beta.cov 6.04e-05
## x0.X1 -5.73e-01
## x0.X2 -5.90e-01
## x0.X3 -7.16e-01
## x0.X4 -6.62e-01
## x0.X5 -5.82e-01
## x0.X6 -6.32e-01
## x0.X7 -5.60e-01
## x0.X8 -4.34e-01
## x0.X9 -7.52e-01
## x0.X10 -6.55e-01
## x0.X11 -6.77e-01
## x0.X12 -5.62e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1618 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1618 iterations.
## Log-likelihood: -138.2772
## AIC: 298.5543 AICc: 300.6833
##
## Estimate
## R.diag 0.210
## Q.alpha.var 0.338
## Q.alpha.cov 0.338
## Q.beta.var 0.097
## Q.beta.cov 0.097
## x0.X1 -0.494
## x0.X2 -0.534
## x0.X3 -0.407
## x0.X4 0.189
## x0.X5 0.444
## x0.X6 0.433
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1992 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1992 iterations.
## Log-likelihood: -359.999
## AIC: 753.998 AICc: 756.3519
##
## Estimate
## R.diag 4.97e-01
## Q.alpha.var 1.01e+00
## Q.alpha.cov 1.01e+00
## Q.beta.var 8.00e-05
## Q.beta.cov 4.74e-05
## x0.X1 -4.03e-01
## x0.X2 -5.21e-01
## x0.X3 -5.00e-01
## x0.X4 -3.95e-01
## x0.X5 -4.40e-01
## x0.X6 -5.11e-01
## x0.X7 -4.63e-01
## x0.X8 -5.31e-01
## x0.X9 -2.76e-01
## x0.X10 -6.60e-01
## x0.X11 -3.44e-01
## x0.X12 -7.99e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.7.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)
3.8 v4: a diag & unequal; b diag & unequal
## empty list for model results
mod_res <- vector("list", length(mpg_ID))
## empty beta forecasts
betas4 <- betas4_se <- NULL
## empty alphas & betas
alphas <- alphas_se <- betas <- betas_se <- NULL
## number of regr params = level + slope
mm <- 2
## CAPMs for each MPG
for(i in 1:length(mpg_ID)) {
## print loop ID
print(paste0("i = ",i,"; MPG = ",mpg_ID[i]))
## get popns
assets_mpg <- assets[grep(mpg_ID[i], rownames(assets)),]
## number of assets
nn <- dim(assets_mpg)[1]
## Process eqn
## B is Identity
BB <- "identity"
## U is col vec of 0's
UU <- "zero"
## Q will be block diagonal
QQ <- matrix(list(0),mm*nn,mm*nn)
## upper L block for alpha
aa <- 1:nn
## equal var-cov
# QQ[aa,aa] <- "alpha.cov"
# diag(QQ[aa,aa]) <- "alpha.var"
## diagonal & unequal
diag(QQ[aa,aa]) <- paste0("alpha",seq(nn))
## lower R block for beta
bb <- aa + nn
## diagonal & unequal
diag(QQ[bb,bb]) <- paste0("beta",seq(nn))
## Observation eqn
## Z is array of regr variables (ie, 1's and market returns)
ZZ <- array(NA, c(nn,mm*nn,TT))
for(t in 1:TT) {
ZZ[,,t] <- cbind(1,market[t]) %x% diag(nn)
}
## A is col vec of 0's
AA <- "zero"
## assume same obs var & cov
##RR <- matrix(list(0),nn,nn)
##RR[aa,aa] <- "obs.cov"
##diag(RR[aa,aa]) <- "obs.var"
RR <- "diagonal and equal"
## Initial states
##x0 <- rbind(matrix("pi.a",nn,1), matrix("pi.b",nn,1))
##V0 <- 100*diag(mm*nn)
##x0 <- rbind(matrix(0,nn,1), matrix(1,nn,1))
##V0 <- matrix(list(0),mm*nn,mm*nn)
##diag(V0[aa,aa]) <- "var.a"
##diag(V0[bb,bb]) <- "var.b"
## starting values for regr parameters
ini_list <- list(x0=matrix(1, mm*nn, 1))
## list of model matrices & vectors
mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR)
##mod_list <- list(B=BB, U=UU, Q=QQ, Z=ZZ, A=AA, R=RR, x0=x0, V0=V0)
## list of control values
con_list <- list(safe=TRUE, allow.degen=FALSE,
conv.test.slope.tol=2, abstol=1e-04, maxit=5e04)
## fit multivariate DLM
mod_res[[i]] <- MARSS(assets_mpg, inits=ini_list, model=mod_list, control=con_list)
## get list of Kalman filter output
kf_out = MARSSkfss(mod_res[[i]])
## ts of forecasted betas
betas4 <- cbind(betas4,t(kf_out$xtt1[(1:nn)+nn,]))
## ts of SE of betas; nxT matrix
betas4_se <- cbind(betas4_se,sqrt(t(apply(kf_out$Vtt1,3,diag)[(1:nn)+nn,])))
## ts of alphas
alphas <- cbind(alphas,t(mod_res[[i]]$states[1:nn,]))
## ts of SE{alphas}
alphas_se <- cbind(alphas_se,t(mod_res[[i]]$states.se[1:nn,]))
## ts of betas
betas <- cbind(betas,t(mod_res[[i]]$states[(1:nn)+nn,]))
## ts of SE{betas}
betas_se <- cbind(betas_se,t(mod_res[[i]]$states.se[(1:nn)+nn,]))
}## [1] "i = 1; MPG = GR"
## Success! abstol and log-log tests passed at 1378 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1378 iterations.
## Log-likelihood: -418.6651
## AIC: 887.3302 AICc: 892.4483
##
## Estimate
## R.diag 8.43e-01
## Q.alpha1 6.52e-05
## Q.alpha2 1.32e-01
## Q.alpha3 2.11e-02
## Q.alpha4 1.01e-01
## Q.alpha5 8.29e-05
## Q.alpha6 5.10e-02
## Q.beta1 2.41e-02
## Q.beta2 1.59e-01
## Q.beta3 1.26e-01
## Q.beta4 2.72e-04
## Q.beta5 4.27e-02
## Q.beta6 7.89e-03
## x0.X1 -6.42e-02
## x0.X2 1.07e+00
## x0.X3 -1.95e-01
## x0.X4 -1.10e-01
## x0.X5 -8.11e-02
## x0.X6 -1.90e-01
## x0.X7 -4.61e-02
## x0.X8 2.25e+00
## x0.X9 -1.07e+00
## x0.X10 -1.11e-01
## x0.X11 -1.62e-01
## x0.X12 -4.89e-02
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 2; MPG = MF"
## Success! abstol and log-log tests passed at 3158 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 3158 iterations.
## Log-likelihood: -458.1922
## AIC: 966.3844 AICc: 971.7341
##
## Estimate
## R.diag 1.45e+00
## Q.alpha1 1.16e-02
## Q.alpha2 4.30e-01
## Q.alpha3 1.04e-04
## Q.alpha4 1.23e-01
## Q.alpha5 4.24e-05
## Q.alpha6 4.23e-05
## Q.beta1 1.07e-04
## Q.beta2 1.01e-04
## Q.beta3 6.94e-05
## Q.beta4 1.11e-04
## Q.beta5 2.75e-02
## Q.beta6 7.91e-02
## x0.X1 -4.03e-01
## x0.X2 -1.31e+00
## x0.X3 -2.40e-01
## x0.X4 -6.99e-01
## x0.X5 -1.67e-01
## x0.X6 -3.34e-01
## x0.X7 -3.17e-01
## x0.X8 -3.57e-01
## x0.X9 -4.62e-01
## x0.X10 -5.31e-01
## x0.X11 -1.90e-01
## x0.X12 -2.84e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 3; MPG = SF"
## Success! abstol and log-log tests passed at 1259 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1259 iterations.
## Log-likelihood: -172.4242
## AIC: 370.8483 AICc: 373.832
##
## Estimate
## R.diag 4.27e-01
## Q.alpha1 1.08e-02
## Q.alpha2 1.59e-03
## Q.alpha3 7.53e-05
## Q.beta1 3.87e-01
## Q.beta2 2.22e-01
## Q.beta3 1.25e-01
## x0.X1 -3.68e-01
## x0.X2 -1.41e-01
## x0.X3 1.64e-01
## x0.X4 6.99e-02
## x0.X5 9.89e-01
## x0.X6 7.69e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.
##
## [1] "i = 4; MPG = SR"
## Success! abstol and log-log tests passed at 1331 iterations.
## Alert: conv.test.slope.tol is 2.
## Test with smaller values (<0.1) to ensure convergence.
##
## MARSS fit is
## Estimation method: kem
## Convergence test: conv.test.slope.tol = 2, abstol = 1e-04
## Estimation converged in 1331 iterations.
## Log-likelihood: -462.4741
## AIC: 974.9483 AICc: 980.107
##
## Estimate
## R.diag 3.92e-01
## Q.alpha1 1.06e+00
## Q.alpha2 5.17e-01
## Q.alpha3 5.60e-01
## Q.alpha4 5.79e-01
## Q.alpha5 5.83e-01
## Q.alpha6 2.42e+00
## Q.beta1 1.53e-04
## Q.beta2 9.44e-05
## Q.beta3 9.76e-05
## Q.beta4 9.17e-05
## Q.beta5 5.42e-01
## Q.beta6 3.28e-04
## x0.X1 -5.04e-01
## x0.X2 -1.05e+00
## x0.X3 -8.87e-01
## x0.X4 -2.71e-01
## x0.X5 -4.24e-02
## x0.X6 6.98e-02
## x0.X7 -4.12e-01
## x0.X8 -5.86e-01
## x0.X9 -3.01e-01
## x0.X10 -4.91e-01
## x0.X11 -1.33e-01
## x0.X12 -8.43e-01
##
## Standard errors have not been calculated.
## Use MARSSparamCIs to compute CIs and bias estimates.## assign names
colnames(alphas) <- colnames(alphas_se) <- names_pop
colnames(betas) <- colnames(betas_se) <- names_pop
colnames(betas4) <- colnames(betas4_se) <- names_pop3.8.1 Plots of results
betas <- betas[,names_mat$ID]
# total num of ts
NN <- dim(betas)[2]
# set up colors; group by MPG
col.pal <- c(rainbow(6, start=0.05, end=0.12),
rainbow(3, start=0.2, end=0.45, v=0.7),
rainbow(6, start=0.8, end=0.95, v=0.85),
rainbow(6, start=0.55, end=0.7))
# num of rows for plot
nr <- 3
# num of cols for plot
nc <- 7
# set graphics params
par(mfcol=c(nr,nc), mai=c(0.1,0.07,0.2,0.07), omi=c(0.6,0.6,0,0))
# get extent of y-limits
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.5)*0.5,
ceiling(max(betas+2*betas_se, na.rm=T)/0.5)*0.5)
yext <- c(floor(min(betas-2*betas_se, na.rm=T)/0.1)*0.1,
ceiling(max(betas+2*betas_se, na.rm=T)/0.1)*0.1)
# draw plots
for(i in 1:NN) {
plot(t_index, betas[,i], type="n",
ylim=yext, xaxt="n", yaxt="n",
xlab="", ylab="", main=names_mat[i,"long"])
rect(par()$usr[1], -1, par()$usr[2], 1, col=gray(0.9), border=gray(0.9))
abline(h=0, lty="dashed")
box()
lines(t_index, betas[,i], col=col.pal[i], lwd=2.2)
lines(t_index, betas[,i]+2*betas_se[,i], col=col.pal[i], lwd=0.7)
lines(t_index, betas[,i]-2*betas_se[,i], col=col.pal[i], lwd=0.7)
#text(xin, yin, names.pop[i], adj=c(0,1))
if(i<=nr | nc==1) {
# axis(side=2, at=seq(ceiling(yext[1]/2)*2,floor(yext[2]/2)*2), las=1)
# axis(side=2, at=seq(yext[1],yext[2])[(seq(yext[1],yext[2])%%1)==0])
axis(side=2, las=1, tick=TRUE)
}
if(i%%nr==0) {
axis(side=1, at=t_index[(t_index %% 10)==0])
}
}
#mtext("Brood year", side=1, outer=TRUE, line=3, cex=1.2)
mtext("Year", side=1, outer=TRUE, line=3, cex=1.2)
mtext(expression(beta[italic(t)]), side=2, outer=TRUE, line=2.2, cex=1.5)